LIGO-P0900002 



ON 
O 
O 

(N 

< 

(N 



o 
cr. 

wo. 



> 
m 

On 

o 

On 
O 



X 
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We assess the statistical errors in estimating the parameters of non-spinning black-hole binaries 
using ground-based gravitational-wave detectors. While past assessments were based on partial 
information provided by only the inspiral and / or ring-down pieces of the coalescence signal, the 
recent progress in analytical and numerical relativity enables us to make more accurate projections 
using "complete" inspiral-merger-ringdown waveforms. We employ the Fisher information-matrix 
formalism to estimate how accurately the source parameters will be measurable using a single 
interferometric detector as well as a network of interferometers. Those estimates are further vetted 
by full-fledged Monte-Carlo simulations. We find that the parameter accuracies of the complete 
waveform are, in general, significantly better than those of just the inspiral waveform in the case of 
binaries with total mass M > 20Mq. In particular, for the case of the Advanced LIGO detector, 
parameter estimation is the most accurate in the M = 100 — 2OOM0 range. For an M = 100A/© 
system, the errors in measuring the total mass and the symmetric mass-ratio are reduced by an 
order of magnitude or more compared to inspiral waveforms. Furthermore, for binaries located at 
a fixed luminosity distance dn, and observed with the Advanced LIGO-Advanced Virgo network, 
the sky-position error is expected to vary widely across the sky: For M = WOMq systems at 
dz, — lGpc, this variation ranges mostly from about a hundredth of a square-degree to about a 
square-degree, with an average value of nearly a tenth of a square-degree. This is more than forty 
times better than the average sky-position accuracy of inspiral waveforms at this mass-range. For 
the mass parameters as well as the sky-position, this improvement in accuracy is due partly to 
the increased signal-to-noise ratio and partly to the information about these parameters harnessed 
through the post-inspiral phases of the waveform. The error in estimating dh is dominated by the 
error in measuring the wave's polarization and is roughly 43% for low-mass (M ~ 20Mq) binaries 
and about 23% for high-mass (M ~ WOMq) binaries located at dt = lGpc. 

PACS numbers: 04.30.Tv,04.30.-w,04.80.Nn,97.60.Lf 



I. INTRODUCTION 

Astrophysical black holes (BHs) are typically classi- 
fied into three groups: stellar-mass BHs (with a mass of 
approximately 3 — 30M Q ), [super] massive BHs (~ 10 4 
- 10 10 M Q ) and intermediate-mass (IM) BHs (~ 30 - 
10 4 M Q ). There is strong observational evidence for the 
existence of both stellar-mass and supermassive BHs. 
The existence of stellar-mass BHs, which are the end 
products of stellar evolution, has been primarily inferred 
from observations of X-ray binaries that allow us to esti- 
mate the mass of the compact object through measure- 
ments of the orbital period and the maximum line-of- 
sight Dopplcr velocity of the companion star [l|. The 
mechanism for producing supermassive BHs is less cer- 
tain but the acceleration of gas disks in the bulges of 
nearly all local massive galaxies point to their existence 
there (2j. Even more convincingly, the observations of 
stellar proper motion in the center of the Milky Way have 
confirmed the presence of a supermassive BH 0] . On the 
other hand, the observational evidence for IMBHs is only 
suggestive. The main hint comes from the observations 
of ultraluminous X-ray sources, combined with the fact 
that several globular clusters show evidence for an excess 
of dark matter in their cores 141 . 
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According to hierarchical galaxy-merger models, [su- 
per] massive BH binaries should form frequently, and 
should be common in the cores of galaxies. There is at 
least one piece of clear evidence for the existence of a su- 
permassive BH binary, namely, the X-ray active binary 
black hole (BBH) at the center of the galaxy NGC 6240, 
which is expected to coalesce in Hubble time [f|. There 
is also growing observational evidence for the existence 
of many other [super] massive BBHs @, 0, H ©]. De- 
spite the lack of any observational evidence for stcllar- 
mass/intermediate-mass BH binaries, different mecha- 
nisms to form these binaries have been proposed in the 
literature (see, for e.g., 0, El El ) • 

Coalescing BH binaries are among the most promis- 
ing sources of gravitational waves (GWs) for the ground- 
based interferometric detectors. What makes them ex- 
tremely interesting is that their gravitational waveforms 
can be accurately modelled and well parametrized by 
combining a variety of analytical and numerical ap- 
proaches to General Relativity. To wit, the gravitational 
waveforms from the inspiral stage of the binary can be 
accurately computed by the post-Newtonian (PN) ap- 
proximation to General Relativity, while those from the 
ring down stage can be computed using BH perturba- 
tion theory. The recent breakthrough [lj, EE EH in nu- 
merical relativity has made it possible to compute accu- 
rate gravitational waveforms from the hitherto unknown 
merger stage as well [II EE EE EE EE El HI HI ■ 

Concomitant with that breakthrough has been the no- 
table progress in GW instrumentation. The Initial LIGO 
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(LIGOI) [22j detectors have completed their first science 
run at design sensitivity. The Virgo detector [23| ran 
concurrently with LIGO for part of that run. Currently, 
both observatories are undergoing commissioning work 
with the target of achieving second-generation sensitivi- 
ties over the next several years, to usher us into the era 
of Advanced LIGO (AdvLIGO) H and Advanced Virgo 
(Adv Virgo). Also, an intermediate, enhanced stage of 
LIGO, called Enhanced LIGO (EnhLIGO), is expected 
to be operational this year. 

In the absence of any observational evidence of stellar- 
mass/intermediate-mass BH binaries, the rate of binary 
coalescence events is estimated by population synthesis 
studies. Plausible rate estimates for stellar-mass BH co- 
alescences detectable by LIGOI / EnhLIGO / AdvLIGO 
detectors range from 7 x 1CV 4 / 7 x 1CV 3 / 2 per year 
to 2 / 20 / 4000 per year with a likely rate estimate of 
around 0.01 / 0.1 / 30 per year HI . For the case of IMBH 
binaries, the plausible rates for LIGOI / AdvLIGO de- 
tectors are 10~ 4 /0.1 per year [l2T| . Similarly, for the 
case of stellar-mass BHs merging with IMBHs (the so- 
called intermediate-mass-ratio inspirals) , plausible event 
rates for LIGOI / AdvLIGO are 10~ 3 / 10 per year [lj] 1 . 
A network of interferometric detectors involving LIGO, 
Virgo, and perhaps others, such as GEO600 [27], and 
TAMA [28|, will be able to extract a host of physical 
parameters of those sources, complementing other de- 
tectors probing their electromagnetic characteristics. 

Indeed, some of the BBH mergers, e.g., triggered 
by the mergers of galaxies/stellar clusters harboring 
supermassive/intermediate-mass BHs, are likely to have 
electromagnetic (EM) counterparts. To associate an EM 
event with a GW signal from such a merger, and vice 
versa, one needs to be able to locate the GW source 
with a high enough accuracy so that the number of 
star clusters or galaxies in the sky-position error box 
is sufficiently small. As argued in Rcf. (29[, even arc- 
minute resolution can make such associations quite fea- 
sible. Whereas the GW observations are expected to 
provide more accurate distance measurements than their 
EM counterpart, the latter will locate the sources in 
the sky with far greater resolution than the former. 
This complementarity was explored in Ref. [30j to ar- 
gue that by combining GW and electromagnetic obser- 
vations it should be possible to constrain the values of 
certain cosmological parameters. In particular, using 
the distance-rcdshift relation from many BBH "standard 
sirens" , such multi-messenger observations can put in- 
teresting constraints on the equation of state of the dark 
energy |3ll . [32j . Supermassive BH binaries are also ex- 
cellent test beds for "strong-field" predictions of General 
Relativity (see, e.g., [H, [13|). Also, GW observations of 
BBH coalescences can be used to test theoretical pre- 
dictions such as the "no-hair" theorem [35[. The effec- 
tiveness of these and other applications depends on the 



accuracy with which we can estimate the parameters of 
the binary, which includes the component masses, dis- 
tance, orientation, and sky location. 

In this work, we study the effect of detector noise 
in limiting the accuracy with which parameters of a 
BBH system can be determined with the present and 
planned earth-based laser interferometers. In the past, 
in the absence of complete coalescence waveforms arising 
from numerical relativity, parameter estimation studies 
were constrained to address this question only for the 
inspiral/ring-down pieces of the sig nal present in the 
band of a detector [Hi, M, H, M, Sfl El, E3, 0, 0, S| . 
Here we extend those studies to estimate how the astro- 
physical quest for characterizing such systems benefits 
from the knowledge of the complete coherent signal, com- 
prising some or all of the inspiral, merger, and ringdown 
pieces, that lies in a detector's observational band. Im- 
provements in the accuracy of BBH parameter measure- 
ments might be expected owing to the increased signal- 
to-noise ratio (SNR) arising from the inclusion of the 
post-inspiral pieces. A second avenue toward parameter 
accuracy improvements can also arise, for some parame- 
ters, from the breaking of some parameter degeneracies 
that the extra information carried by the GW phasing 
of those pieces might offer. We employ the phenomeno- 
logical inspiral-merger-ringdown waveforms proposed in 
Refs. [H, Hi, S3 to explore these possibilities. 2 The 
systematic errors that might arise in observations using 
these "com plet e" 3 BBH coalescence templates are stud- 
ied in Ref. [54[ . 

To estimate the parameter errors, we adopt a two- 
pronged approach. One of these is of obtaining the 
Fisher information matrix and then inverting it to de- 
rive the parameter error variance- covariance matrix 55 1. 
The elements of this matrix are then used to obtain the 
lower bound on the parameter estimator errors [56l . . 
This approach is employed here, in spite of its known 
limitations [37l . l38l . [Bq , since it has been studied exten- 
sively in the community and allows for a fair comparison 
of our results with those given in the literature. How- 
ever, since by its very design, this bound may not be re- 
spected for signals with a low SNR (as first demonstrated 
by Refs. [13, [38|), we also assess estimator errors through 
Monte Carlo studies. For the parameter ranges consid- 
ered here, the latter approach corroborates the findings 
of the former, with a few notable exceptions arising from 
parameter space boundaries, where the Monte Carlo es- 
timates reflect better the results expected from real-data 
searches. 

In addition to addressing the primary question on how 
large the parameter errors are, we also study their be- 
havior across the BBH parameter space. We study how 
the various estimator errors scale with the mass param- 
eters themselves. How much improvement do the com- 
plete waveforms impart to the determination of the sky- 
position of BBHs in multi-detector searches? How does 



It should be noted that these assessments take into account only 
the inspiral stage (for the case of stellar-mass and intermediate- 
mass-ratio binaries) or ring-down stage (for the case of IMBH 
binaries) of the binary coalescence. The event rates are likely 
to be higher for a search using inspiral-merger- ring down tem- 
plates. See, for example, Fig. 14 of [2t| for a comparison of the 
sensitivities of searches employing different templates. 



2 A similar study using the effective-one-body-numerical-relativity 
waveforms [4a. |49|. [50 . 51, 52] is being pursued as well |53l . 

3 Throughout this paper, we refer to the waveforms modelling 
all the three (inspiral, merger and ringdown) stages of BBH 
coalescence as "complete" waveforms. 
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the sky-position accuracy change with the BBH mass pa- 
rameters? A summary of our results is as follows: First, 
we find that the parameter-estimation accuracies using 
the complete waveforms are, in general, significantly bet- 
ter than those using only their inspiral phases in the case 
of BBHs with a total mass M = (mi + m 2 ) > 20M Q , 
where are the component masses, at least for mass- 
ratios between 0.25 and unity. The observed trend sug- 
gests that this improvement can be expected for some- 
what lower mass-ratios as well. Second, for BBHs at 
a fixed effective distance and M > WMq whereas the 
fractional errors in the two mass parameters, M and 
r\ = mim^/M 2 , scale mostly monotonically with M for 
the inspiral-only waveforms, they do not display that 
property for the complete waveforms. In the latter case, 
they instead exhibit a distinct minimum, whose location 
is determined by M, 77, and the detector's noise power 
spectral-density (PSD). Third, owing to the use of com- 
plete vis a vis inspiral-only waveforms the sky-position 
accuracy improves by factors of many. We also show that 
for the complete waveforms alone, the sky-position ac- 
curacy mostly degrades with increasing total-mass when 
the SNR is kept fixed. This is primarily caused by a simi- 
lar degradation in the estimation accuracy of the signal's 
times of arrival at the different detectors in a network. 
This deterioration in accuracy, while not monotonic in 
M at finer scales, is broadly so at large scales, and is 
caused by the reduction in the number of in-band wave 
cycles. 

More specifically, for Advanced LIGO, the estima- 
tion of the total mass, the symmetric mass-ratio, and 
the effective distance d g is the most accurate in the 
M = 100 - 2OOM range. (For other detectors, that 
mass range is somewhat different since it is partly de- 
termined by their noise PSDs.) For such systems, the 
reduction of errors in parameter estimates is by an ordcr- 
of-magnitude or more due to the inclusion of the post- 
inspiral phases. The improvement is mainly due to the 
expected increase in SNR arising from the inclusion of 
those phases. This expectation, which is based on the 
assumed Gaussianity and stationarity of detector noise, 
must be tempered by the observation that the amount 
of increase in SNR can be less in real data. 

We also observe that for a fixed SNR, the inclusion 
of the post-inspiral phases improves the accuracy of M 
and 77 for a wide range of masses much more (by several 
times) than that of the chirp mass M c . This is due to the 
fact that the inclusion of those phases helps in breaking 
the degeneracy between those two parameters (M and 
77) known to exist in the inspiral waveform. 

For a fixed SNR, the estimation of the luminosity dis- 
tance for low-mass systems shows negligible change by 
the inclusion of the post-inspiral phases. This is due to 
its strong covariance with the polarization and the or- 
bital inclination angles of the binary, which is mostly 
unaltered by the inclusion of the post-inspiral phases. 
Also, for a fixed SNR, the luminosity distance estimate 
deteriorates with increasing M, for reasons discussed be- 
low. On the other hand, for a fixed luminosity distance, 
the error in its estimate initially improves with increas- 
ing M, due to the increase in SNR, before degrading 
eventually owing to the decreasing number of in-band 
wave cycles. 

Before moving on, we wish to point out some limi- 



tations of the present work. First, this study considers 
only the dominant harmonic of non-spinning BBH wave- 
forms. Astrophysical BHs are expected to have spin, and 
including spin effects can change the estimation of dif- 
ferent BBH parameters (59|- Whereas on the one hand 
previous calculations have shown that the parameter- 
estimation accuracies generally deteriorate upon the in- 
clusion of spin-orbit and spin-spin couplings [36l . |59| . 
on the other hand the inclusion of spin-induced preces- 
sion in the waveform model can improve the parame- 
ter estimation [f5(| [6l|. Also, it has been noted in var- 
ious studies that including the higher harmonics can 
significantly increase the parameter-estimation accura- 
cies [H Hi, il |H, [H Hi]. So, while the results pre- 
sented in this paper may not be too far from the real- 
istic case, we stress that a rigorous statement on the 
parameter-estimation accuracies should consider these 
effects as well. Moreover, neglecting spins and higher 
harmonics in the waveform models can result in signif- 
icant amount of systematic errors in estimating various 
parameters. These systematic errors are out of the scope 
of this paper. A preliminary investigation of this is pre- 
sented in Ref. |54| . 

This paper is organized as follows: Sec. |TT] briefly in- 
troduces the main aspects of the search for binary black 
holes. In particular, Sec. Ill Al reviews the phenomenolog- 
ical inspiral- merger-ring down waveform templates pro- 
posed in Refs. [H, Hfl while Sec. Ill Til provides a 
brief introduction towards the statistical theory of pa- 
rameter estimation. In Sec. IIHI we present the results 
of our calculations in the case of a search using a sin- 
gle interferometric detector. This section discusses the 
results from the analytical calculations using the Fisher- 
matrix formalism as well as the numerical Monte-Carlo 
simulations. Results from the calculations in the case 
of a network of detectors are discussed in Sec. IIVI while 
Sec. |V] summarizes the main results and provides a dis- 
cussion of the possible astrophysical implications of this 
work. 



II. GRAVITATIONAL WAVE OBSERVATIONS 
OF NON-SPINNING BINARY BLACK HOLES 

In General Relativity, the gravitational-wave strain at 
any point in space can be expanded in terms of its two 
linear polarization components h + (t) and h x (t) or the 
two related circular polarization components, 

W) = h+(t)-ihx(t)=A(t)j*® (2.1) 

and its complex conjugate, with ip(t) and A(t) denoting 
the wave's phase and amplitude. Generally, the GW 
emitted by a coalescing binary has multiple harmonics. 
In this work, we limit our study to only the dominant 
harmonic's contributions to tp(t) and A(t). Then the GW 
strain h(t) in a detector is the linear combination of the 
two polarization components, h(t) = F + h + (t)+F x h x (t), 
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FIG. 1: Noise amplitude spectrum (\/Sh(f)) of different de- 
tectors considered in this paper. 



with the detector's antenna-pattern functions given as: 

F + (6,<p,4>) = --(1 + cos 2 0) cos 2(f) cos2V> 

— cos 9 sin 20 sin 2ip, 
1 



F x (6,<j),4>) = -(1 + cos 2 6) cos 20 sin2*0 
— cosO sin 20 cos 2ip. 



(2.2) 



Above, 9 and <j) arc the polar and azimuthal angles speci- 
fying the location of the source in the sky in the detector 
frame and ip is the polarization angle. 

The two polarization components of the BBH signals 
are sinusoids with varying amplitude and frequency, and 
have phases 7r/2 radians apart relative to each other. 
Consequently, their GW signal in a detector can be writ- 
ten as: 



h(t) =CA(t) cos[cp(t) +cp ], 



(2.3) 



where the amplitude coefficient C and phase ipo can be 
assumed to be constant for signals lasting for a dura- 
tion (up to several minutes) much shorter than Earth's 
rotational time-scale: 



C 



1 



(1 + cos 2 i) 2 Fl + 4 cos 2 lF 



tan 



2F X cos i 
F+(l + cos 2 i) 



(2.4) 



Above, i is angle of inclination of the orbit to the line of 
sight. 



A. Detecting non-spinning binary black holes 

The GW signal's phase (f(t) and amplitude A(t) are 
functions of the physical parameters of the binary, such 
as the component masses and the spins. Detecting a 
signal requires analyzing intcrferometric data, which are 



noisy. Defining a search strategy, therefore, necessitates 
the modelling of this noise, which we take here to be 
zero-mean Gaussian and stationary: 

Mtj = 0, (2.5) 
n*(/)n(/') = \s h {f)5{f-f) , (2.6) 

with the over-bar denoting the ensemble average and the 
tilde denoting the Fourier transform, 



n(f) 



n{t) e 



-2m ft 



dt. 



(2.7) 



Above, Sh(f) is the Fourier transform of the auto- 
covariance of the detector noise and is termed as its (one- 
sided) power spectral-density. We also assume the noise 
to be additive. This implies that when a signal is present 
in the data x(t), then 



x(t) = h(t) + n(t) . 



(2.8) 



The noise covariance Eq. (|2 . 6[) introduces the following 
inner-product in the function space of signals: 



(a, b) 



a*(f)kf) 



(2.9) 



where a(f) and b(f) are the Fourier transforms of a(t) 
and b(t), respectively. 

Under the above assumptions about the characteris- 
tics of detector noise, the Neyman-Pearson criterion [55[ 
leads to an optimal search statistic, which when max- 
imized over the amplitude coefficient C, is the cross- 
correlation of the data with a normalized template, 



p = (h, x) , 



(2.10) 



where the normalized template is h(f) = h(f)/ \J (h, h). 
In a "blind" search in detector data, where none of the 
binary's parameters are known a priori, the search for 
a GW signal requires maximizing p over a "bank" of 
templates (see, for e.g., (67|) corresponding to different 
values of those physical parameters. Apart from the 
physical parameters, the waveform also depends on the 
(unknown) initial phase (po and the time of arrival to. 
Maximization over the initial phase ipo is effected by us- 
ing two orthogonal templates for each combination of the 
physical parameters [68j . and the maximization over t 
is attained efficiently with the help of the Fast Fourier 
Transform (FFT) algorithms [69j . 

Since the cross correlation between the data and the 
template can be most efficiently computed in the Fourier 
domain by using the FFT, waveform templates in the 
Fourier domain are computationally cheaper. Refer- 
ence [26| proposed a family of analytical Fourier domain 
templates for BBH waveforms of the form: 



(2.11) 



where the effective amplitude and phase are expressed 
as: 
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A eS (f) 

*eff(/) 



^ /merg ^ (///merg) /0 if /mcrg < / < /ring 
«>£(/) /ring, Cr) if /ring < / < /cut, 
1 7 

27r/i + + - V(x fc t? 2 + VkV + Zk) (irMf)^- 5 ^ 3 . 

77 * — ' 

' fc=0 



(///merg)" 7/6 if / < /„ 



-2/3 



(2.12) 



r 



In the above expressions, 



1 



£(/, /ring, 0-) = 



2^7 (/ - /ring) 2 + ^ 2 /4 



(2.13) where A 



is a Lorentzian function that has a width cr, and that 

The normal- 

2/3 " 

, is chosen so as 



is centered around the frequency / r ; ng 



ization constant, w = ^ f J"" g ^ 



to make A c s(f) continuous across the "transition" fre- 
quency /ring- The parameter / mcr g is the frequency at 
which the power-law changes from /~ 7 / 6 to /~ 2 / 3 . The 
effective distance to the binary is denoted by c? e ff, which 
is related to the luminosity distance <1l by d e ff = cLl/C. 
The phcnomenological parameters /merg, /ring, cr and / cu t 
are given in terms of the total mass M and symmetric 
mass-ratio 77 of the binary as 



a i] 2 + b Q i] + c Q , 



7rM/ mcrg 

TT-W/ring 

7tM(J = a 2 7] 2 + 62 + C 2 

7rM/ cut = a 3 7? 2 + 63 77 + c 3 



(2.14) 



The coefficients a,j,bj,Cj, j = 0...3 and Xk,VkjZk, k 
0, 2, 3, 4, 6, 7 are tabulated in Table I of Ref. El. 



B. Measuring binary black hole parameters 

To evaluate how effective the detectors will be in estab- 
lishing the field of GW astronomy, especially, with the 
second-generation Earth-based interferometers sched- 
uled to come online around 2014, one needs to foremost 
assess how accurately they can measure the astrophysical 
properties of compact object binaries. That quest will be 
limited, on the one hand, by the accuracy with which the 
search templates can model actual gravitational wave- 
forms, and, on the other hand, by the inherent statisti- 
cal noise in the measurement process. The former issue 
is one of systematics, which will be discussed elsewhere 
(see, e.g., Ref. [54]). Here, we discuss the latter issue in 
more detail. 

To determine how large the noise-limited errors can 
be in the measured values of the signal parameter, we 
take those values to be the maximum likelihood estima- 
tors (MLEs). The discussion in the preceding section 
shows that a total of nine parameters characterize the 
non-spinning BBH coalescence signals considered here. 
They are the total mass M, the symmetric mass-ratio 
77, the sky-position angles (9, c6), the binary's orienta- 
tion angles (ip,i), the luminosity distance cLl, the initial 
(or some reference) phase tpo, and the time of arrival 
(or some reference time) t$. For computing the error 



estimates, we map them onto the components of the 
parameter vector, # = {In .4, to, <po, InM, In 77, 9, 4>, ip, 

M 5/6 fEn 



Owing to noise, their MLEs, 



d cff7 r 3 /2 Y 24- 

i9, will expectedly fluctuate about the true values, i.e., 
= i9 + <5i?, where 5i9 a is the random error in estimat- 
ing the parameter 7? a . The magnitude of these fluctua- 
tions can be quantified by the elements of the variance- 
covariance matrix, 7 ab = ch9 a Sd b (5f|. 

A relation between the j ab and the signal is available 
through the Cramer-Rao inequality, which dictates that 



7 



> 



(2.15) 



where T is the Fisher information matrix: 
[d a h I ('d),d b h I (€f)) 



ab 



N 

E 

7=1 



(-0 



I=1 J a h\J) 

where I is the detector index and d a denotes taking par- 
tial derivative with respect to the parameter i? a . There- 
fore, A7? a = ( 6$ a 5$ a ) V2 = Taa 2 gives the lower 
bound on the root-mean-square (rms) error in estimat- 
ing i) a . The two are equal in the limit of large SNR (see, 
e.g., S3). 

The errors in the sky-position angles will be presented 
in terms of the error in the measurement of the sky- 
position solid angle, defined as: 



AO = 2ir\J(Acos9A(f>) 2 



5 cos 9 



(2.17) 



Each parameter-error, Ad a , falls off inversely with SNR. 
Since the solid angle is two dimensional, its error falls off 
quadratically with SNR [H, [zO, [HI ■ 



III. PARAMETER ESTIMATION: 
SINGLE-DETECTOR SEARCH 

A. Analytical calculation using Fisher information 
matrix 

In this section, we use the Fisher information-matrix 
formalism to estimate the errors in measuring the param- 
eters of coalescing BBHs with a single GW interferome- 
ter. We present results for three generations of ground- 
based detectors, namely, Initial LIGO, Enhanced LIGO 
and Advanced LIGO. The one-sided noise PSD of the 
Initial LIGO detector is given in terms of a dimension- 
less frequency x = f / fo by [72|, [73[ 




FIG. 2: Errors in estimating the total mass M (top left), symmetric mass ratio r\ (top middle), chirp mass M c (top right), 
time of arrival to (bottom left) and effective distance d e « (bottom right) in the case of Advanced LIGO noise spectrum, plotted 
against the total mass of the binary. The errors of M, r/, M c and d e ff are in percentage and the errors of to are in seconds. 
The value of the symmetric mass-ratio rj is shown in legends. The solid lines correspond to a search using complete BBH 
templates and the dashed lines correspond to a search using 3.5PN-accurate post-Newtonian templates in the SPA, truncated 
at the Schwarzschild ISCO. The binary is placed optimally oriented at an effective distance of 1 Gpc. 



S h {f{x)) = 9 x 10~ 46 [(4.49a;)- 56 + 0.16a;- 4 - 52 + 0.52 + 0.32a; 2 ] , 
where /o = 150 Hz; while the same for Enhanced LIGO reads [74| : 

S h (f(x)) = 1.5 x 10- 46 [l.33 x iO" 27 e - 5 - 5 ( ln:E ) 2 a;- 52 - 6 + 0.16 a;- 4 - 2 + 0.52 + 0.3 a; 21 

where f = 178 Hz. For Advanced LIGO [l2], 

S h (f{x)) = 10- 49 

where /o = 215 Hz, and, for Advanced Virgo [7 
S h (/(*)) = lO" 47 



x" 4 - 14 - 5a;~ 2 + 111 



1 - a; 2 + a; 4 /2 
1 + a; 2 /2 



2.67 X 10- 7 X;- 5 - 6 + 0.59 e< ln ^ ["3.2-1.08 ln(,)-0.13(ln X -A.l + gg e -0.73(ln X f ^ 



(3-1) 



(3.2) 



(3.3) 



(3.4) 



where / = 720 Hz. The calculations presented in 
this section were performed using the Initial LIGO, En- 
hanced LIGO and Advanced LIGO noise PSDs, while 
the calculations presented in Sec. IIVI consider a three- 
detector network consisting of Advanced LIGO and Ad- 
vanced Virgo. 

The parameters that can be estimated through single- 
detector observations are {A(<1l), to, (fo, M, 77}. To 
be precise, one can measure only the Doppler-shiftcd 
masses, unless there are additional experiments for de- 
termining the Doppler-shift [36l | and, therefore, allow the 



estimation of the true masses. Doppler shifting can arise 
due to the motion of the detector relative to the source 
or the cosmological expansion. In measurements with 
multiple-detectors, as discussed below, it is possible to 
measure the source distance and sky-position as well. 
There too, the distance observed is actually the Doppler- 
shifted distance. 

The Fisher matrix elements in the {.A, io, fo, M, 77} 
space are computed from the derivatives of the wave- 
forms described by Eqs. ([2~TT|) - ([2~14|) : 
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FIG. 3: Same as in Fig. [2] except that the binary is placed at an effective distance of 100 Mpc and the noise PSD corresponds 
to that of Initial LIGO. 
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FIG. 4: Same as in Fig. [2] except that the binary is placed at an effective distance of 100 Mpc and the noise PSD corresponds 
to that of Enhanced LIGO. 



Tab = {d a h(f),d b h(f) 



, , d a A cS (/) d b A cS (/) + A 2 cfi (f) dAff (/) gftgeffCf) 



(3.5) 



r 



where the low- frequency cutoff, /i OW i is chosen to be 10 
Hz for Advanced LIGO, and 40 Hz for Enhanced and 
Initial LIGOs. The upper-frequency cuttoff. / cut is given 



by Eq. (ICT) . 

The rms errors in parameters M, rj and to are com- 
puted by inverting the Fisher matrix elements as dis- 
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FIG. 5: The overlap function (i.e., the ambiguity function maximized over to and <fio) between waveforms constructed at 
different points in the parameter space. The horizontal axis reports the total mass M of the binary while the vertical axis 
reports its symmetric mass-ratio n. Each panel shows the overlap of different waveforms with one "target waveform". The 
total mass of the target waveform is chosen to be M — 20Mq, IOOMq, 200Mq, 400Mq, respectively, for the four columns 
starting from the left. The symmetric mass-ratio of the target waveforms is chosen to be n = 0.25, 0.222, 0.16 in the top, 
middle and bottom rows, respectively. 



cussed in Sec. Ill B[ The error in estimating the chirp 
mass M c and the effective distance d e f{ are obtained by 
propagating the errors in M, r\ and A in the following 
way: 
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(3.7) 



where A-d a denotes the rms error in estimating $ a ob- 
tained from T a b, and C a b is the correlation coefficient 
between parameters $° and r d b . 

Errors in the estimates of the parameters M, r), M c , to 
and d e s in the case of AdvLIGO detector arc plotted 
against the total mass M in Fig.[2j These errors are com- 
puted assuming that the binary is placed at an effective 
distance of 1 Gpc. Also plotted in the figures are the 
same error-bounds computed from the 3.5PN accurate 



restricted PN waveforms in the stationary phase approxi- 
mation (SPA), truncated at the Schwarzschild innermost 
stable circular orbit (ISCO). It can be seen that, over 
a significant range of the total mass, the error-bounds 
in the complete templates are largely better than those 
in the PN inspiral waveforms. For binaries with M = 
IOOA/q and r] = 0.25, the error-bounds in various param- 
eters using the complete [PN] templates are AM/M ~ 
0.34 [5.38]%, Arj/r] ~ 0.84 [12.98] %, AM C /M C ~ 
0.35 [2.47]%, A£ ^ 0.46 [15.51] ms and Ad e $/d eS ~ 
1.36 [5.24] %. The errors in estimating the same param- 
eters using Initial LIGO and Enhanced LIGO detectors 
are plotted in Figs. [3]and|3J 

The rate of variation in the errors in different regions 
of the parameter space can be understood by studying 
the overlap function, which is the ambiguity function 
maximized over to and ifo [zE] ■ Figure [5] plots the con- 
tours of the overlap between waveforms generated at dif- 
ferent points in the (M, 77) space. Notice the change in 
the shape and orientation of the ambiguity ellipses, es- 
pecially, as the total mass of the binary is varied. While, 
to a very good approximation, the chirp mass contin- 
ues to remain as one of the eigen-coordinate [77| in the 
case of the low-mass (with M < 20M Q ) binary inspi- 
ral (PN) waveforms, this is no longer true for the com- 
plete waveforms of higher mass systems. This is be- 
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FIG. 6: Errors in estimating the total mass M (top left), symmetric mass ratio rj (top middle), chirp mass M c (top right), 
time of arrival to (bottom left) and effective distance d e g (bottom right) in the case of Advanced LIGO noise spectrum, plotted 
against the total mass of the binary. The errors of M, rj, M c and d e ff are in percentage and the errors of to is in seconds. 
Symmetric mass-ratio 77 is shown in legends. The solid lines correspond to a search using complete BBH templates and the 
dashed lines correspond to 3.5PN-accurate post-Newtonian templates in the SPA, truncated at the Schwarzschild ISCO. The 
errors correspond to a fixed SNR of 10. 




FIG. 7: Same as Fig. [6] except that the noise PSD corresponds to that of Initial LIGO. 



cause the latter waveforms have more information about 
the component masses than just the chirp mass. The 
eigen-dircctions change dramatically with increasing to- 
tal mass. It can be seen that the error trends reported 
in Fig. [2] closely follow the shape of these ambiguity el- 



lipses. This also means that while placing templates in 
the inspiral-merger-ring down searches, we will have to 
consider these changes in the orientation of the ambigu- 
ity ellipses. This will be studied in a future work. 
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One common problem encountered in the estimation 
of errors using Fisher information matrix is the follow- 
ing: In some cases (especially in the case of large number 
of parameters), the Fisher matrix becomes badly condi- 
tioned, thereby, decreasing the fidelity of the error co- 
variance matrix derived by inverting it. This problem 
can often be obviated by intelligently choosing the pa- 
rameters and by projecting out certain dimensions in the 
Fisher matrix (e.g., to and ipo)- We have verified our re- 
sults by comparing the errors computed using the full 
Fisher matrix with those computed using the projected 
matrix. In our calculations, they turned out to be the 
same to the extent discernible in the figures and tables 
presented here. 

It may be noted that for a fiducial signal limited only 
to the inspiral phase of the binary, i.e., for / < f merg , 
the parameter A is uncorrelated with the other signal pa- 
rameters, and hence one has Ti a = Si a p 2 , which renders 
the Fisher matrix in the block-diagonal form. However, 
for the complete signal, with the merger and the ring- 
down pieces included, the correlation of A with the other 
parameters becomes non zero, and the Fisher matrix is 
no longer block-diagonal with respect to this parameter. 
This implies that the complete waveforms provide more 
information about A and, hence, about the effective dis- 
tance (i e ff. 

Figures [7] and [5] show the error estimates corre- 
sponding to a fixed (single-detector) SNR of 10 in the 
case of Advanced LIGO, Initial LIGO and Enhanced 
LIGO noise spectra, respectively. It is interesting to note 
that the parameter estimation using the complete wave- 
forms is still much better than that using only the inspi- 
ral waveform even though, in order to produce the same 
SNR using inspiral templates, the effective distance to 
the binary has to be often much smaller. The reason for 
this can be understood through an analogy with param- 
eter estimation with multiple detectors: Since the inspi- 



ral phase, on the one hand, and the merger-ringdown 
phases, on the other hand, occupy two contiguous and, 
essentially, non-overlapping frequency-bands, the detec- 
tion of a complete signal is equivalent to a coherent de- 
tection of these two pieces of the waveforms by two coin- 
cident, co-aligned detectors with sensitivities limited to 
the two contiguous bands, respectively. The two phases, 
however, are modulated by the two mass parameters in 
complementary ways, in the sense that the Fisher sub- 
matrices in the two-dimensional mass-space for these two 
fiducial detectors grow more linearly independent of each 
other, the larger the total mass gets, even while the total 
coherent SNR of this fiducial detector pair is held con- 
stant. This linear independence causes the estimation of 
two mass parameters to improve. Contrastingly, since 
the merger-ringdown pieces add very little information 
about a system's chirp-mass, the improvement in its ac- 
curacy arising from using complete waveforms is much 
less even for high mass systems. 

Figure[5]plots the SNR produced at different detectors 
by equal-mass binaries located at a fixed distance, as a 
function of the total mass of the binary. 



B. Monte-Carlo simulations 

The limitations of the Fisher-matrix formalism are 
well known [3?], [H, |58| . The parameter-error bounds 
provided by it are trustworthy in the limit of high SNR 
and for parameters on which the signal has linear depen- 
dence. In the case of low SNRs the error bounds com- 
puted using the Fisher matrix formalism can be largely 
different from the "actual" errors. Also, the Fisher ma- 
trix does not recognize the boundaries of the parameter 
space (such as the restriction r\ < 0.25). Neither does it 
account for practical restrictions such as the finite sam- 
pling of the data. In order to explore these limits of 
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3.80 
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1.11 
0.51 
0.18 



4.40 
3.13 
1.54 
0.69 
0.29 



10.0 
6.02 
3.01 
1.50 
0.60 



18.8 
11.3 
5.64 
2.82 
1.13 



5.59 
3.57 
1.94 
0.95 
0.31 



6.93 
4.96 
2.93 
1.48 
0.52 
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0.32 
0.13 
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(4.64) 
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(4.22) 
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(12.0) 


10 


0.19 
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(2.62) 


11.9 


(5.56) 


0.17 


(0.35) 2.60 


(3.20) 


11.4 
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20 


0.10 
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2.60 
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7.25 
5.26 
3.01 
1.45 
0.53 



8.50 
5.10 
2.55 
1.28 
0.51 



20.2 
12.1 
6.06 
3.03 
1.21 



6.58 
4.49 
2.73 
1.33 
0.46 



10.2 
9.58 
6.43 
3.25 
1.16 



1.66 
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0.50 
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0.10 



3.06 
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4.58 
2.58 
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0.61 
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9.98 
6.64 
3.44 
1.69 
0.67 



5.11 
3.07 
1.53 
0.77 
0.31 



23.6 
14.2 
7.07 
3.54 
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Ad c ff/rfeff X 100 



6 


16.7 


17.6 


24.6 


16.7 


17.5 


24.3 


16.7 


17.4 


23.2 


10 


10.0 


10.5 


14.8 


10.0 


10.5 


14.6 


10.0 


10.4 


13.9 


20 


5.00 


5.27 


7.39 


5.00 


5.26 


7.30 


5.00 


5.21 


6.95 


40 


2.50 


2.63 


3.70 


2.50 


2.63 


3.65 


2.50 


2.61 


3.48 


100 


1.00 


1.05 


1.48 


1.00 


1.05 


1.46 


1.00 


1.04 


1.39 


Ato(ms) 


6 


0.37 


5.90 


15.8 


0.39 


7.22 


20.0 


0.42 


11.7 


36.3 


10 


0.22 


3.54 


9.47 


0.23 


4.33 


12.0 


0.25 


7.04 


21.8 


20 


0.11 


1.77 


4.73 


0.12 


2.17 


6.00 


0.13 


3.52 


10.9 


40 


0.06 


0.88 


2.37 


0.06 


1.08 


3.00 


0.06 


1.76 


5.45 


100 


0.02 


0.35 


0.95 


0.02 


0.43 


1.20 


0.03 


0.70 


2.18 



TABLE I: Errors in estimating different prameters of the binary, as a function of the binary parameters and the SNR, computed 
using the Fisher matrix formalism. Same errors computed from the Monte-Carlo simulations are shown in brackets. 



the Fisher formalism, we performed Monte-Carlo sim- 
ulations, whereby maximum- likelihood detections were 
made of simulated signals added to multiple statistically 
independent realizations of simulated colored, Gaussian 
noise. The aim of this frequentist study was to obtain the 
spread in the maximum-likelihood estimates of the pa- 
rameters and compare them with Fisher-matrix calcula- 
tions. It is worth clarifying that there is another interest- 
ing question one can pose in the context of parameter es- 
timation, namely, "Given a specific signal and a particu- 
lar noise realization, what are the posterior distributions 
of the parameter estimates." This is a question from 
Bayesian statistics that can be answered using Markov- 
Chain Monte-Carlo (MCMC) simulatio ns, as explored 
for inspiral-only waveforms in Refs. [ll, [TiJ HH . We 
do not answer that question here. 

In this section we present results from the frequen- 
tist Monte-Carlo simulation studies. These studies 
largely corroborate the Fisher matrix calculations in the 
parameter-space regions where the latter is expected to 
be trustworthy. The simulations also allow us to com- 
pute error-bounds in the parameter-space regions where 



the Fisher matrix formalism can be unreliable (such as 
for rj ~ 0.25). We caution the reader that this is not 
meant to be an exhaustive comparison between Fisher- 
matrix calculations and Monte-Carlo simulations. A 
detailed comparison of Fisher matrix formalism with 
Monte-Carlo simulations in the case of 3.5PN inspiral 
signals can be found in the recent work Ref. [HJ . 

Colored Gaussian noise with one-sided PSD Sh(f) is 
generated in the frequency domain. If x k and y k de- 
note the real and imaginary parts of the discrete Fourier 
transform of the noise at the frequency bin fc, these are 
generated by 

Xk = \J~Sh~ k x k /2 , y k = y r S^y k /2, (3.8) 

where x k and y k are random variables drawn from a 
Gaussian distribution of zero mean and unit variance, 
and Sh k denotes the discrete version of Sh{f )- Frequency 
domain signal described by Eq. (|2.11| ) is added to the 
noise. The data is filtered through a matched filter em- 
ploying templates described by Eq. (|2.1ip . The likelihood 
is maximized over to and ipo as described in Sec. [Ill The 
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FIG. 9: The curves labelled "IniLIGO" and "EnhLIGO" re- 
port the SNR produced by binaries located at an effective 
distance of lOOMpc at Initial LIGO and Enhanced LIGO, re- 
spectively, as a function of the total mass. The curves labelled 
"AdLIGO" and "AdVirgo" report the same produced by bi- 
naries located at lGpc at Advanced LIGO and Advanced 
Virgo. The solid lines correspond to complete waveforms and 
the dashed lines correspond to PN waveforms. 
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FIG. 10: Scatter plot of parameters estimated from 10 4 
Monte-Carlo simulations. The horizontal axis reports the 
total mass and the vertical axis reports the symmetric mass 
ratio. The left panel correspond to the injection with pa- 
rameters M = 20Af© and r\ = 0.2222, and the right panel 
correspond the injection with parameters M = 200Mq and 
rj = 0.2222. The injections correspond to an SNR of 20. Also 
overlaid in the left panel is a cartoon of the four different ini- 
tial simplexes chosen for the maximization algorithm. The 
true values of the parameters is marked by a cross. Note that 
the eigen directions are different in the two plots. 



maximization over the physical parameters (M and rf) 
is best performed by filtering the data, using a template 
bank finely spaced in the parameter space. But, in order 
to attain sufficiently good accuracy (say, 1%), a large 
number of simulations needs to be performed. Thus, 
computing error-bounds from a good volume of the pa- 
rameter space is computationally expensive in a template 
bank search. So, in this paper, the maximization over 
the physical parameters is performed with the aid of the 
computationally cheaper Nelder-Mead downhill simplex 
algorithm 

We emphasize that this search may not be as accu- 
rate as the template bank search. One reason for the 



inaccuracy is that, in this method, we do not "sample" 
the parameter space finely enough, and hence the "real 
maximum" can very well be missed. This is especially 
the case when the function that we want to maximise 
(likelihood in this case) contains many secondary max- 
ima. Indeed, it is well known that the the likelihood can 
have many secondary maxima arising due to global cor- 
relations in the parameter space. We bypass this issue 
by starting the maximisation algorithm around the "ac- 
tual" peak of the function. Hence, the error distributions 
that we obtain are only indicative of the spread of the 
MLM estimates around the primary maxima. Unlike in 
the case of MCMC simulations, this does not provide a 
complete picture of the posterior distribution of the pa- 
rameters. Nevertheless, this is a worthwhile tool as an 
independent verification of the Fisher matrix calculation, 
enabling us to "scan" ' a good volume of the parameter 
space using Monte-Carlo simulations 4 . 

Neldcr-Mead's algorithm is a multidimensional min- 
imisation/maximisation algorithm. In order to maxi- 
mize the required function, we need to specify an initial 
"simplex" of n + 1 dimensions where n is the dimension- 
ality of the parameter space. Since the dimensionality 
of our parameter space is 2, the simplex in our case is 
a triangle. It is important for the good convergence of 
the maximization that the initial simplex "catch" the 
orientation of the ambiguity ellipses in our parameter 
space, which often depends strongly on the parameters 
themselves. Thus, we start the maximization by speci- 
fying four different initial simplexes, whose vertices have 
equal (coordinate) distance from the "true" value of the 
parameters. The four triangles arc oriented in different 
directions in the parameter space. We choose the param- 
eters corresponding to the best among the maximized 
likelihoods as the parameters of the injection. Figure [TU] 
shows a scatter plot of the parameters estimated from 
10 4 simulations. Also overlaid in the left plot is a car- 
toon of the initial simplexes chosen. The reader may note 
the difference in the eigen-directions in the two plots. 

We found that the following points need to be taken 
care of while performing this kind of simulations: (i) 
Since the frequency-domain templates are abruptly cut 
off at the frequency / cu t, we need to make sure that 
the edges arising from this do not corrupt our numerical 
calculations. This means that, for high mass systems 
(M > 200-Mq) we cannot perform the simulations with 
very high SNR (p > 100), because the cutoff frequency is 
at the "sweet spot" of the detector, (ii) Sufficiently small 
tolerance level for the maximization algorithm in order 
to ensure that the "true" maximum is never missed, (iii) 
Orthonormality of the search templates, as emphasized 
by Ref. 

The frequency distributions of the estimated param- 
eters M, rj and M c are shown in Fig. [11] The injection 
corresponds to the parameter values M = 20M© and 
77 = 0.16 and an SNR of 20. Also plotted in the fig- 
ures are the expected distributions computed using the 



4 In our simulations, a few hundred trials were sufficient for the 
Nelder-Mead's algorithm to converge to the fiducial maximum. 
By contrast, a template bank search requires tens of thousands 
of templates, in general. 
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FIG. 12: The data points show the errors computed from Monte-Carlo simulations in the case of Advanced LIGO noise PSD. 
The horizontal axis reports the total mass while the legends report the symmetric mass ratio. The errors are computed for a 
fixed SNR of 20. The dashed lines correspond to the same errors computed using Fisher matrix formalism. 



Fisher matrix formalism. All the results are computed 
using the AdvLIGO noise PSD. It can be seen that the 
two calculations agree very well. Figure [T^] shows the er- 
rors computed using the Monte-Carlo simulations plot- 
ted against the total mass of the binary for three differ- 
ent values of 77. The simulations are performed with an 
SNR of 20. Also shown are the error-bounds computed 
using the Fisher matrix formalism. In the case of mass 
ratios 77 = 0.2222 and 77 = 0.16, the simulations agree 
well with the Fisher matrix calculations. But the simu- 
lations disagree with the Fisher calculations for the case 
of rj = 0.25. This is expected because the Fisher matrix 
does not recognize the physical restriction that 77 can 
only take values less than, or equal to 0.25. The Fisher 
matrix calculation assumes that the errors in estimat- 
ing the parameters are Gaussian distributions centered 
around 77 = 0.25, while the Monte-Carlo simulations en- 
force the restriction 77 < 0.25. As a result the error 
bounds estimated by the Monte-Carlo simulations will 
be less than that estimated by the Fisher matrix. 

Fisher matrix calculations assume that the errors de- 
crease inversely proportional to the SNR. But this ap- 
proximation is not valid at low SNRs. So we have per- 
formed Monte-Carlo simulations with various SNRs in 
order to study the SNR dependence of the errors. Fig- 
ure [TJJ] plots the errors estimated from the simulations 
against the SNR of the injections. The top, middle and 
bottom panels in the figure correspond to mass ratios 
77 = 0.25,0.2222 and 0.16, respectively. The different 
markers correspond to the Monte-Carlo simulations and 
the dashed lines correspond to the Fisher matrix calcu- 



lations. It can be seen that, barring the case of rj = 0.25, 
the simulations agree very well with the Fisher calcula- 
tions in the limit of high SNRs (p > 10). Because of 
the 77-boundary effects, the errors computed from the 
r/ = 0.25 simulations are less than those computed from 
the Fisher calculations. For small SNRs (p < 10), the 
simulation errors start to deviate from the Fisher calcu- 
lations. There are two reasons for this: (i) at low SNRs, 
as observed by many others (see, for e.g., Ref. (38|) the 
Fisher matrix largely underestimates the errors. This 
is the dominating effect in the case of M = 2OM bi- 
naries at low SNRs in Fig. [13] (ii) at low SNRs, since 
the size of the ambiguity ellipses are increased, they are 
cut by the 77 = 0.25 boundary, which is neglected by 
the Fisher calculations. Hence the Fisher matrix over 
estimate the errors. This is the dominating effect in the 
case of M = 200M Q binaries at low SNRs. It is the inter- 
play between these two competing effects that causes the 
discrepancy between the simulations and Fisher calcula- 
tions. In summary, the results from the Monte-Carlo 
simulations, albeit the limitations of the maximization 
algorithm used, should be more reliable than the Fisher 
calculations. 



Table [I] tabulates the errors in the case of Advanced 
LIGO noise PSD, computed using both Fisher matrix 
and Monte-Carlo simulations. 




FIG. 13: Errors computed from Monte-Carlo simulations (crosses, dots and diamonds) plotted against SNR. The horizontal 
axes report the SNR of the injections and the legends report the total mass in units of Mq. The top, middle and bottom 
panels correspond to mass ratios n = 0.25,0.2222 and 0.16, respectively. The error-bounds expected from the Fisher matrix 
calculation are indicated by dashed lines. 



IV. PARAMETER ESTIMATION: 
MULTI-DETECTOR SEARCH 



With a sufficiently large number of geometrically in- 
dependent and well-separated interferometric detectors 
it is possible to measure all nine of the BBH parameters 
of an adequately strong source [H, [H| . To assess how 
accurately such a measurement can be made with the 
AdvLIGO-AdvVirgo network, one can begin by comput- 
ing the Fisher matrix in the nine-dimensional parameter 
space, and then invert it to obtain the error variancc- 
covariance matrix. We take the network to comprise 
three interferometers, with one each at Hanford (WA), 
USA, Livingston (LA) USA, and Cascina, Italy. The 
LIGO detectors in Hanford and Livingston are assumed 
to be having the AdvLIGO noise PSDs given in Eq. (|3.3p 
and the Virgo detector in Cascina is assumed to be hav- 
ing the Adv Virgo noise PSD given in Eq. (|3.4[) . 

When interpreting the astrophysical implications of 
these parameter errors, it is important to remember that 
it is only when the signal is linear in the parameters or 
the SNR is large that the maximum-likelihood estima- 
tor is unbiased and the error deduced from the Fisher 
matrix achieves the Cramer-Rao bound 



55(. To aid this 



conformity, we map four of the six extrinsic signal pa- 
rameters (i.e., parameters that depend on the observers 
location in time and space), viz., (A, tp, h Vo); into new 



parameters, a k , with k =1,...,4, such that the signal in 
Eq. (|2.3[) at any given detector has a linear dependence 
on them: 



h(t) =^a fe h fc (t) 



(4.1) 



fe=i 



where the hfc(£)'s are completely independent of those 
four extrinsic parameters. (The two remaining extrinsic 
parameters are the sky-position angles.) To deduce their 
dependencies as well as the forms of the a fc 's we begin by 
noting that the antenna-pattern functions can be treated 
as the components of a vector that are related to two sky- 
position dependent functions, u(9, <p) and v(8, <p) (4ll.[8q. 
through a two-dimensional rotation by 2tp: 



Fy, 



cos 2tp sin 2-0 
— sin 2ip cos 2ip 




(4.2) 



With this well-known observation, one finds 

hi(£) oc u(6,<p) cos[(p(t)] , 

h a (t) cx v(e,<f>)cos[<p(t)}, 

h 3 (i) oc u(6, 4>) sm[ip{t)] , 

h 4 (t) cx v(6,4>)sm[p(t)], 



(4.3) 



where the proportionality factor is a dimensionlcss 
(mass-dependent) function of time. 
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FIG. 14: The network SNR of a signal, corresponding to 
the complete waveform, from an equal-mass binary with 
M — IOOA/0 located at dh = lGpc, plotted as a function 
of its sky-position. The network here is the three detector 
AdvLIGO-AdvVirgo network, such that the two 4km-arm- 
length LIGO detectors in Hanford and Livingston have Ad- 
vLIGO noise PSDs and the Virgo detector in Cascina has 
AdvVirgo noise PSD. Above, 9 and (f> are the polar and az- 
imuthal angles specifying the location of the source in the 
sky in the geographic coordinate system. 



The new parameters are themselves defined as 
i a 1 a 3 \ _ y{i) 



M a 



where j/(t) 



O vo ■ X ■ 02-i/i 



(4.4) 



'l + cos 2 /,) + 4 cos 2 /, 



1/2 



O a is the two- 



dimensional orthonormal rotation matrix for angle a and 



T = f (1 + cos 2 t) /y(i) 

\ 2cosi/y(t) 



(4.5) 



The Fisher matrix is then computed on the space 
(M, r], 8, cj>, to, a 1 , a 2 , a 3 , a 4 ). The errors in the a fc 's are 
obtained by inverting that matrix. By using error- 
propagation equations obtained from Eq. (|4.4[) . we are 
able to deduce error estimates for all four extrinsic pa- 
rameters. 

In this paper, however, we present the error esti- 
mates for, perhaps, the most astrophysically interesting 
of those, namely, the luminosity distance. To obtain it, 
first notice that 



tv{Mj M a ) = \\a\\ 



4 



(4.6) 



where tr is the trace, and ||a|| 2 = ^fe=i i ak ) 2 ■ This 



yields 



d(d L 



dy 
V 



dllall 



(4.7) 



which can then be used to deduce the rms error, Adi/d^ 
by accounting for the covariance between y and a . Fi- 
nally, we choose a flat prior in y(i), such that whenever 



its estimate is negative or greater than its maximum pos- 
sible value (of four) the prior is set to zero. The distance 
errors plotted below are for such a prior. 

The error variance-covariance matrix described above 
can also be used to derive the error estimates for the 
other astrophysically interesting quantity, namely, the 
sky-position. Here again, to further keep our assess- 
ment robust, we first reduce the dimensionality of the 
Fisher matrix to five by projecting out the four above- 
mentioned extrinsic parameters. This helps in lower- 
ing the condition number of the Fisher matrix across 
the parameter space. We do so by taking a cue from 
Refs. [H|, HH, where it was shown that the network like- 
lihood ratio of compact binary inspiral signals can be 
maximized analytically over those four extrinsic param- 
eters. Moreover, just as for the signal in a single detec- 
tor, it is possible to speed up the search in to by using 
the FFT [85|. Thus, the only parameters that need to 
be searched numerically through the help of a template 
bank [8(| are the following four parameters: (M, 77, 6, <f>). 

The resulting Fisher matrix is well-behaved every- 
where in the five-dimensional sub-space except on a set 
of points of measure zero, where the detectors in the net- 
work cease to be geometrically independent. Its inverse 
yields the error estimates for the two mass parameters 
and the sky-position. A sky-map of the network SNR 
is presented in Fig. 1141 while the sky-maps of the errors 
in the source luminosity-distance and the sky-position 
are given in Fig. [TS] for equal-mass BBH sources with 
M=100M o and located at d L = lGpc. 

Figure [TBI shows the all-sky distribution of the errors 
in estimating the solid angle f2. The left plots show 
the probability density and the right plots show the cu- 
mulative distribution. We assume that the sources are 
distributed uniformly across the sky. Top panels cor- 
respond to a binary with M = 20M Q and 77 = 0.25, 
and the bottom panels to a binary with M = 100M Q 
and r) = 0.25. In each plot the thick (red) traces corre- 
spond to the errors estimated using the complete wave- 
forms while the thin (black) traces correspond to those 
estimated using restricted 3.5PN waveforms in the SPA 
truncated at Schwarzschild ISCO. All the errors are com- 
puted for a network SNR of 10 for the respective wave- 
forms. The error-estimates are obtained by averaging 
over the angles (ijj, t). These plots show that in the case 
of an M = 2OM and r\ = 0.25 binary, assuming that 
the sources are distributed uniformly across the sky, the 
sky-position of 70% [10%] of the sources can be estimated 
with an accuracy better than 1 [0.1] square degree. Us- 
ing PN templates, the sky-location of only 29% [6%] of 
the sources can be estimated with an accuracy better 
than 1 [0.1] square degree. For the M = 100M Q binary, 
the sky-position of 90% [18%] of the sources can be es- 
timated with an accuracy of 1 [0.1] square degree using 
complete waveforms, while only 15% [4%] of the sources 
can be resolved with the same accuracy using inspiral 
waveforms. It should be noted that in that figure we 
have normalized the errors for SNR fixed to 10. For real 
systems additional improvement might be seen from the 
use of the complete waveforms provided their inclusion of 
merger and ringdown phases actually improves the SNR 
of those signals. This is indeed the case for high-mass 
systems (M > 20M Q ). For an equal-mass binary with 
M = 20 [100] M Q , the improvement in the SNR by the 
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FIG. 15: The left plot shows the sky-position error log 10 [ASl (in square-degrees)] and the right plot shows the fractional error 
in the luminosity distance log 10 [Ad_L/d_L (in %)] as functions of the sky-position of a BBH source. The source studied here is 
the same equal-mass binary considered in Fig. 1141 and, 9 and c/> are the polar and azimuthal angles specifying the location of 
the source in the sky in the geographic coordinate system. Note how the effect of the varying network sensitivity, as seen in the 
SNR plot in Fig. 1141 is imprinted in the two error plots. Additionally, the error plots display a full "sine-wave" pattern, which 
comprises a set of sky-positions for which the geometric independence of the LIGO- Virgo detectors is the weakest. Extraction 
of the signal's polarization is affected the most at these locations. That in turn hurts the distance measurement accuracy. 
The same locations do not necessarily hurt the determination of the sky-position, which is mostly driven by the measurement 
accuracy of the times-of-arrival of the signal at the three sites. 
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FIG. 16: All-sky distribution of errors in estimating the solid 
angle in the case of AdvLIGO-AdvVirgo network. The 
left plots show the probability density and the right plots 
show the cumulative distribution. The top panels correspond 
to an equal-mass binary with M = 20Mq, and the bottom 
panels to one with M = WOMq. In each plot the thick (red) 
traces correspond to the errors estimated using the complete 
waveforms while the thin (black) traces correspond to those 
estimated using restricted 3.5PN waveforms. All the errors 
are computed for a network SNR of 10 for the respective 
waveforms. 
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FIG. 17: All-sky distribution of errors in estimating c?l in the 
case of AdvLIGO-Adv Virgo network. The left plots show the 
probability density and the right plots show the cumulative 
distribution. The top panels correspond to an equal-mass 
binary with M — 20Mq, and the bottom panels to one with 
M = IOOMq . In each plot the thick (red) traces correspond 
to the errors estimated using the complete waveforms while 
the thin (black) traces correspond to those estimated using 
restricted 3.5PN waveforms. All the errors are computed for 
a network SNR of 10 for the respective waveforms. 



inclusion of merger and ringdown is 9% [300%] , in sta- 
tionary, Gaussian noise. (See the discussion of Fig. [TBI 
below.) 

Figure [T7] shows the distribution of the errors in esti- 



mating the luminosity distance to two different types 
of equal-mass binary systems, both producing a network 
SNR of 10 in the AdvLIGO-AdvVirgo network. The 
top panels correspond to a binary with M = 20M Q and 
the bottom panel to a binary with M = IOOMq. As 
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FIG. 18: Same as Fig. [16] except that the binary is now placed 
at a fixed luminosity distance of lGpc. Notice the strong sim- 
ilarity between the plots in the top panel above and those in 
the top panel of Fig. \W\ This is because in the plots of the 
top panel above the average SNR is relatively close to 10. 
The plots in the bottom rows of the two figures are more dis- 
parate: The average SNR above is several [few] times better 
than the fixed SNR in Fig. [16] for the complete [inspiral-only] 
waveforms. 
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FIG. 19: Same as Fig. [17] except that the binary is placed 
at a fixed luminosity distance of lGpc. By comparing the 
above figure with Fig. 1171 it is manifest that nearly all the 
improvement in the luminosity-distance measurement accu- 
racy, when including the post-inspiral phases, arises due to 
the increased SNR. 



in Fig. [Tni the thick (red) traces correspond the com- 
plete waveforms while the thin (black) traces correspond 
to the PN waveforms. These plots suggest that for an 
SNR of 10 the luminosity distance to around 10%[50%] 
of the sources can be estimated with an accuracy of bet- 
ter than 38% [53%] in the case of low-mass systems. They 
also reveal that, for a fixed value of the network SNR, 



the error estimates using inspiral and complete wave- 
forms are almost identical. This is not surprising be- 
cause for low-mass systems the signal is dominated by 
the inspiral phase. In the case of high-mass systems 
with an SNR of 10, the luminosity distance to around 
10% [50%] of the sources can be estimated with an accu- 
racy of 60% [100%]. These errors are worse than those 
for the PN waveform 5 primarily because the covarianccs 
between the initial phase and (i/j, l) are stronger in the 
case of complete waveforms. This property of the com- 
plete waveforms mitigates the estimation accuracy of l, 
which, in turn, affects the estimation of cLl- 



Figures [TBI and [191 show the errors in estimating and 
d,L in the case of binaries distributed uniformly across the 
sky but located at a luminosity distance of lGpc. These 
errors also are averaged over ip and i. These plots show 
that in the case of an equal-mass binary with M = 20M Q 
the sky-position of around 10% [50%] of the sources can 
be estimated with a resolution of 0.07[0.5] square de- 
gree or better. In the case of a M = 100M Q binary, 
10% [50%] of the sources can be estimated with a resolu- 
tion of 0.01 [0.1] square degrees. These plots in Fig. [18] 
also show that the coherent addition of the merger and 
ringdown phases brings about remarkable improvement 
(i.e., by several times for most sky-positions) in the es- 
timation of Q. 



Figure [111] shows that the luminosity distance of 
10% [50%] of the M = 20M Q BBH sources can be es- 
timated with 32% [47%] accuracy or better and that of 
10%[50%] of the M = 100 Af© binaries can be estimated 
with an accuracy of 13% [20%] or better. While compar- 
ing Figs. ITTl and flUl it may help to track the mean errors 
listed in Table [TTJ Studying these plots and numbers 
reveals some interesting aspects of these signals. First, 
for the PN waveforms the distance error improves only 
slightly in going from an SNR of 10 to a source distance 
of lGpc. This is easily explained by the fact that the 
sky-averaged SNR of these systems at di=lGpc is only 
slightly greater than 10. Second, the distance error re- 
duces a little for complete waveforms vis a vis inspiral 
ones at lGpc. This is mainly due to the increased SNR of 
the former. Third, the error for the complete waveforms 
for the M = lOOAf© system at lGpc is still the small- 
est of all the cases studied here because its sky-averaged 
SNR is sufficiently large; indeed, it is large enough to 
even compensate for the increased covariance between ipo 
and (ip, l) arising from the merger and ringdown phases, 
as discussed above. 



5 Note that, in order to get the same SNR in the case of PN 
waveforms, the binary must be placed at a much closer distance. 
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p= 10 


cLl — lGpc 


M/Mq 


AQ Ad L /d L 


AO. Ad L /d L 


20 


0.78 (2.2) 55.7% (55.3%) 


0.70 (2.1) 43.2% (46.8%) 


100 


0.55(8.9) 111% (63.1%) 


0.13(5.9) 23.0% (39.8%) 



TABLE II: Sky-averaged errors in estimating Q and d^ using 
complete BBH waveforms in the case of AdvLIGO-Adv Virgo 
network. The left column tabulates the errors corresponding 
to a fixed value p — 10 for the network SNR, while the right 
column tabulates the errors corresponding to a fixed value 
= lGpc of the luminosity distance. Errors computed 
using PN templates are shown in parentheses. The Q errors 
are given in square degrees and the fractional dz, errors are 
given in percentage. 



Finally, we compare our results with a couple of past 
studies in the form of Refs. [H, |4l|. First, both these 
early studies used the same noise PSD for both LIGO 
and Virgo detectors. Second, their noise PSD was dif- 
ferent from both the AdvLIGO and the AdvVirgo noise 
PSDs used here; it made their detectors more sensitive 
(by a factor of a few in amplitude) in the band below 
70Hz and somewhat less sensitive at higher frequencies 
than the AdvLIGO PSD used here. Third, they consid- 
ered only inspiral signals from binary neutron stars with 
a component mass of IAM@, and distributed them uni- 
formly across a spatial volume. Fourth, in Ref. [4l[ the 
authors culled every source that gave a distance error of 
greater than 100% or that had an SNR of less than 8.5. 
In our study, where all sources were kept at a fixed dis- 
tance of lGpc, none of them were culled. Also, whereas 
all our sources with M = 100 Mq have an SNR greater 
than about 25, those with M = 20 M & have the small- 
est SNR equal to 6. These differences make it difficult to 
compare these different studies. It is, however, possible 
to make some limited comparisons. Specifically, Fig. 15 
in Ref. [4l[ suggests that the fractional errors in the esti- 
mated source distances all tend to be greater than 100% 
as their source distance approaches lGpc. Figure 14 of 
Ref. [36| depicts a similar trend. This appears to be 
consistent with our numbers. 



V. SUMMARY 

In this paper, we studied the statistical errors in es- 
timating the parameters of non-spinning BH binaries 
using ground-based GW observatories. Our study was 
restricted to the leading harmonic of the GW polariza- 
tions of such sources; but employing waveforms mod- 
elling the inspiral, merger and ring-down stages of the 
binary coalescence. We obtain results both for singlc- 
and multi-detector searches. The single-detector prob- 
lem was investigated in the context of two generations 
of ground-based detectors, namely, Initial LIGO and Ad- 
vanced LIGO, as well as Enhanced LIGO, with interme- 
diate sensitivity. On the other hand, the multi-detector 
problem was investigated in the context of the Advanced 
LIGO- Advanced Virgo network. For these calculations, 
we adopted a two-pronged approach: We first analyti- 
cally computed the error bounds using the Fisher-matrix 
formalism. We then pointed out the limitations of this 



approach and improved upon those calculations by full- 
fledged Monte-Carlo simulations. 

To summarize, we find that with an Advanced LIGO 
detector the total mass of an equal-mass binary with 
M = 20Mq[100Mq] located at 1 Gpc can be estimated 
with an accuracy of ~ 0.67[~ 0.34]%, while its sym- 
metric mass ratio can be estimated with an accuracy 
of ~ 1.26[~ 0.84]%. The effective distance can be es- 
timated with an accuracy of ~ 4.87[~ 1.36]% and the 
time-of-arrival can be placed within ~ 0.11 0.46] ms. 
We considered binaries with three different mass ratios 
(t? = 0.25,0.2222,0.16) in the range 1OM < M < 
450M Q for these calculations. These results predict for 
a significantly more accurate astrophysical characteriza- 
tion than what has been presented in the past literature 
(which use the post-Newtonian waveforms intended to 
model only the inspiral stage of the binary). To wit, the 
error-bounds for total mass, computed using the com- 
plete waveforms is better than those computed using the 
inspiral-only waveforms by a factor of ~ 1.4 [~ 16] for an 
equal-mass binary with total mass 2OAf0[lOOM©]. The 
error-bounds on the symmetric mass ratio is improved 
by a factor of ~ 1.4 [~ 15], those on the time-of-arrival 
is improved by a factor of ~ 7 [~ 34] and those on the 
effective distance is improved by a factor of ~ 1.1 [~ 4] 
by the inclusion of the merger and ringdown stages. 

In the case of a network consisting of two Advanced 
LIGO detectors and one Advanced Virgo detector, we 
found that the luminosity distance to an equal-mass bi- 
nary with M = 20M Q at lGpc can be estimated with 
a sky- and orientation-averaged accuracy of 43.2% and 
the sky location can be estimated with a mean accu- 
racy of 0.7 square degrees. For a similar binary, but 
with M = IOOMq, the respective mean accuracies are 
23% and 0.13 square degrees. For low-mass binaries, 
with (M ~ 20M Q ), the improvement in the sky-position 
accuracy due to the inclusion of merger and ringdown 
is about a factor of 3, while for high-mass binaries 
(M r-j 100M Q ), that improvement is by a factor of 45. 
The inclusion of the same two phases betters the distance 
estimates by a few (for low mass systems) to several (for 
high mass systems) percent. In short, the sky resolution 
is greatly improved by the inclusion of merger and ring- 
down, while the improvement in the estimation of the 
luminosity distance arises largely from the extra SNR 
contributed by the merger and ringdown. 

In the case of the AdvLIGO-AdvVirgo detector net- 
work, the parameter-estimation accuracy peaks for bi- 
naries with M ~ 100M Q . Although the observational 
evidence for BHs in this mass range is only suggestive, 
there is growing consensus in the astronomy community 
that IMBHs could exist in dense stellar clusters. The ex- 
istence of this class of black holes could explain a number 
of observations, such as the ultraluminous X-ray sources 
and the excess dark matter concentration in globular 
clusters. 

Several authors have considered the scenario of the co- 
alescence of IMBHs and have come up with coalescence- 
rate predictions [ill, EH ■ Particularly interesting is the 
case of the merger of two stellar clusters each hosting 
an IMBH considered in Ref. |ll|. Since this is expected 



to be a strong source of GW signal with a possible EM 
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counterpart 6 , it is a worthwhile question to ask what 
kind of constraints can be put on the values of cosmolog- 
ical parameters by combining GW and EM observations 
of such sources [3 Of . The improved parameter estimation 
might help to tighten these constraints. This is being in- 
vestigated in an ongoing work [13] • 
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